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INTRODUCTION 

The  differential  equations  of  notion  of  the  Moon  about  its  center 
of  mass  have  been  known  since  the  tine  of  Euler  (Moutsoulas,  1971), 
but  due  to  the  complexity  of  the  driving  terms  (torques  from  Earth,  Sun, 
etc.)*  they  can  only  be  solved  approximately.  Formulations  in  which 
the  orbital  motions  are  approximated  by  functions  containing  terms 
secular  and  periodic  in  time  are  amenable  to  exact  solution,  and  the 
theories  of  motion  derived  are  called  analytic  theories.  Modem 
analytic  theories,  such  as  those  of  Eckhardt  (1970)  and  Migus  (1976) 
are  invaluable  for  their  concise  description  of  different  modes  of 
physical  libration,  and  for  the  possible  detection  of  free  lunar  libra- 
tions,  but  are  constrained  in  accuracy  by  their  dependence  upon  (relatively 
inaccurate)  analytic  orbit  theories. 

In  order  to  obtain  the  improved  accuracy  necessary  for 
the  interpretation  of  modern  lunar  laser  ranging  and  very-long- 
baseline  interferometry  observations,  Williams  (1975)  and  others 
at  the  Jet  Propulsion  Laboratory  developed  a series  of  lunar 
libration  models,  the  latest  being  called  LLB-5,  based  upon 
direct  numerical  integration  of  the  equations  of  motion  for  the 
libration.  The  LL3-5  integration  is  done  in  non-inertial  coordinates, 
so  that  the  equations  of  motion  contain  "inertial-force"  terms. 
Also,  no  variational  equations  are  integrated  in  LL3-5  and  the 
partial  derivatives  needed  for  adjustment  of  the  parameters  of 
the  model  are  obtained  by  finite-differencing. 
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Our  lunar  rotation  is  based  upon  the  direct  numerical 
integration  of  Euler's  differential  equations  for  the  rotation 
of  the  moon,  with  the  rotation  being  described  by  Euler  angles 
referenced  to  an  inertial  coordinate  system.  The  variational 
equations  for  the  partial  derivatives  of  the  Euler  angles  and 
their  time-derivatives  with  respect  to  all  of  the  parameters 
in  the  equations  of  motion,  and  with  respect  to  the  initial 
conditions  of  motion,  are  integrated  in  parallel.  The  detailed 
formulation  of  our  model,  and  results  of  comparing  it  with 
LLB-5,  are  described  in  the  following  sections. 

EQUATIONS  OF  NOTION 


The  state  of  rotation  of  a rigid  body  with  respect  to  an  arbitrary 
coordinate  system  can  be  expressed  by  six  quantities:  three  angles 

defining  a rotation  of  axes,  and  their  temporal  derivatives.  We  choose 
to  describe  the  lunar  rotation  in  terms  of  Euler  angles  and  their 
derivatives  as  defined  by  an  inertial  coordinate  system  referred  to 
the  mean  equinox  andequator  of.  1950.0.  The  axis  is  perpendicular 
to  the  mean  equator  of  1950.0  and  points  northward,  the  axis  is  the 
intersection  of  the  mean  equator  and  ecliptic  of  1950.0  and  points 
towards  the  constellation  Aries,  and  the  ^ axis  completes  the  right- 
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handed  system.  Letting  the  x^  (i  = l,3)  axes  be  a right-handed  system 
coinciding  with  the  Moon's  principal  axes  or  inertia,  where  the  moment 
about  the  x^  axis  is  least  and  the  moment  about  x^  greatest,  the  Euler 
angles  are  defined  as  indicated  in  Figure  1. 

Euler's  equations  governing  the  motion  of  a rigid  body  about  its 
center  of  mass  are 

Atl>1  + (C  - B)u)2u),  = 

B<L2  - (C  - A>1u3  = T2 

Cu>3  + (B  - A) = T,  (1) 

In  this  representation  the  (i  = l,3)  are  the  components  of  the  lunar 
angular  velocity  vector  3 along  the  three  principal  (body- fixed]  lunar 
axes,  x^;  A,  B,  and  C are  the  moments  of  inertia  about  the  principal 
axes;  and  the  T\  are  the  components  of  the  total  torque  vector  about  the 
corresponding  axes. 

Since  we  desire  an  ultimate  accuracy  of  .01"  for  our  model,  we  should 
examine  all  torques  which  might  result  in  displacements  of  this  magnitude. 
The  largest  libration  term  is  due  to  the  regression  of  the  node  of  the 
lunar  equator  on  the  ecliptic  and  is  of  roughly  S000"  amplitude.  Thus, 
so  long  as  we  are  not  driving  near  a resonant  frequency,  we  should  expect 
that  all  torques  that  are  no  more  than  2 x 10  ^ times  the  dominant  term 
(central-body  term  for  the  Earth)  should  be  negligible.  The  far-field 
torque  exerted  on  the  Moon  by  an  external  body  is  directly  proportional  to 
the  body's  mass,  and  inversely  proportional  to  the  cube  of  its  distance. 


fMre.i  n 
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A calculation  of  M/r  . . for  the  bodies  of  the  solar  system  demonstrates 

minimum 

that  the  torque  on  the  Moon  may  be  approximated  as  the  sum  of  the  Earth  and 
Sun  induced  torques  to  an  accuracy  of  1 part  in  10^,  with  the  solar  contribu- 
tion only  1/200  that  of  the  Earth.  The  torque  due  to  the  oblateness  of  the 
Earth  is  not  included  in  the  model  described  in  this  paper,  but  we 
intend  to  incorporate  it  in  the  near  future.  Introducing  the  lunar 


moment  of  inertia  ratios 

a = 

C-B  a C-A  . B-A 

, S * -y  , and  Y = -y  , we 

then  have 

= 

-°“2“3  * (TS1  * V/A 

"2 

= 

8“l“o  * * T02)/B 

“3 

= 

-*V2  * (T93  ♦ XS3)/C 

(2) 

where  and  represent  the  components  of  the  torque  due  to  Earth 

and  Sun,  respectively. 

We  desire  the  equations  of  motion  in  terms  of  the  inertially- 
referenced  Euler  angles.  The  relation  between  the  body-fixed  angular 
rates  and  the  inertial  Euler  angles  is  well  known  (Goldstein,  19S0)  : 

' U1 

= 

9-cos  <p  * sin  <P  sin  6 

U2 

— 

• • 

- 8 sin  b + b cos  b sin  9 

W3 

= 

b cos  9 + b 

(3) 
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Differentiating  the  equations  (3)  with  respect  to  time  results  in  a 
system  of  equations  linear  in  $,  8,  and  $,  the  solution  of  which  forms 
our  Euler  angle  equations  of  motion: 

$ = esc  8 (co1  sin  <t>  + u>2  cos  p)  + 8<i  esc  9 - p6  cot  9 = 

V • • • • 

9 = cos  <p  - u2  sin  p - ip<p  sin  9 = F2 

$ s u3  - cos  8 + ii8sin9  = (4) 

• 

The  ok  are  found  in  terms  of  Euler  angles  by  equations  (2)  and  (3) 
and  it  only  remains  to  evaluate  the  body-axis  torques  due  to  Earth  and 
Sun.  The  potential  due  to  the  lunar  gravitational  field  may  be  approx- 
imated as  a finite  expansion  in  spherical  harmonics  as  follows  (Kaula,  1966) : 


3M  q t 20 

uO,L.  - 1 - 2 J (|)n  P (sin  L) 

r ( n=2  n r n 


where  G = Gaussian  gravitational  constant 
M<j  = mass  of  the  moon 
a = lunar  radius 

r ■ distance  from  the  lunar  center  of  mass 

L * selenocentric  latitude 

X = selenocentric  East  longitude 

Pn  are  the  Legendre  polynomials 

Pnn  are  the  associated  Legendre  functions 

J , C , and  S are  dimensionless  coefficients 
n nm  nm 


In  terms  of  the  rectangular  coordinates  x^^  it  is  easily  shown  that: 


r J 

sin  L = — 


cos  A = 


sin  A = 


In  our  model  the  Earth  and  Sun  are  treated  as  point-masses  of  mass 
(b  s ©,  0)  and  therefore  experience  the  force 

?b  = ’W*  W 

where  V^U  is  the  gradient  of  U evaluated  at  the  position  of  body  b. 
The  resultant  torque  vector  acting  on  the  moon  is 


Tb  = rx(-F) 

= M^xV^)  (8) 

with  r the  vector  from  the  center  of  mass  of  the  Moon  to  the  body.  If 
the  selenocentric  components  of  r are  (Xj , x2,  x3)  then  the  selenocentric 
torque  components  for  body  b are: 
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= M (x  1!L  . V UL 

^ 1 8x2  x2  3Xj 


In  the  evaluation  of  3U/3x^  we  note  that  we  may  neglect  the  central  force 
2 

(l/r  ) term  since  it  produces  no  net  torque  about  the  center  of  mass.  In 
fact  any  term  in  3U/3x^  of  the  form  g*x^,  where  g is  any  function  not 
explicitly  containing  the  index  k,  will  cancel  in  the  torque  cross 
product.  We  obtain  for  either  Earth  or  Sun: 


^ = GHa  nil®  ^ L)  T^T  * GM‘  J2  7 


3 sin  L 


*{[C  cosmX  + S sin  mX]  P'  (sin  L) 

nm  nm  J nm  v 1 3x 


+ m[-C  smmX  + S cosmX]  P (sin  L)  1^—  } 
nm  nm  J nm  v J 3x, 

k 

+ terms  which  do  not  contribute  to  the  torque 
From  equations  (6a)  and  (6b)  we  derive 


3 sin  L 53k  x3xk 


(10) 


(11a) 


where  is  the  Kronecker  delta  defined  by  <$„  = 1 for  i = j and  0 
otherwise.  If  we  differentiate  (6d)  with  respect  to  x^  we  get 


3A  , 2 2„ -3/2  , 2 , „ „ 

9xk  *1  + ^*1  ^2k  ~ *1*2  5lk 


12  °lk' 
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which  by  (6c)  gives: 


X2  Slk 
2 


(lib) 


The  three  second-order  equations  of  motion  are  now  rewritten  as  a 
set  of  six  first-order  differential  equations  for  ease  of  integration 
method.  Defining 


yl 

zz 

hi 

•'t 

X 

$ 

y2 

= 

e 

y5  = 

• 

0 

y3 

r 

y6  E 

$ 

we  obtain : 

• 

yl 

= 

y4 

X • 

-r*. 

ii 

F1 

y2 

= 

y5 

• 

y5  = 

F2 

y3 

s 

y6 

• 

y6  = 

F3 

(12) 

where  the  F. 

l 

are 

given  by  equations 

■ (4). 

The 

six  initial  conditions  of 

the  state  vector  y are  the  three  Euler  angles  and  their  rates  at  some 
initial  epoch. 


VARIATIONAL  EQUATIONS 

Our  primary  reason  for  modelling  the  lunar  physical  librations  is 


to  enable  us  to  process  accurate  data  and  consequently  arrive  at  improved 
estimates  for  the  parameters  affecting  libration.  Often  this  is  accom- 
plished through  the  iterative  use  of  a linear  least-squares  estimator; 
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thus,  in  addition  to  the  model  for  motion,  it  is  desirable  to  supply 
partial  first-derivatives  of  the  state  vector  with  respect  to  the  libra- 
tion  parameters  at  all  times.  Referring  back  to  the  defining  equations 
(4)  for  the  driving  terms  Fk  we  note  that  in  general  the 

Fk  5 Fk^*  P»  y) 

are  explicit  functions  of  time,  the  set  of  adjustable  parameters  p,  and 
the  Eulerian  state  vector  y.  Ke  differentiate  the  equations  of  motion  (12) 
with  respect  to  any  specific  time-independent  parameter  p^  and  interchange 
the  order  of  differentiation  to  obtain  the  variational  equations: 


_d_/^k\  _ 3yk+3 

dt  y3p ~ 3p^  k = 1,3 


(13) 


_d_ 

dt 


t,p  i 


with  p t i denoting  the  parameter  set  p exclusive  of  pi-  The  initial 
conditions  for  these  equations  are  all  zero  except  when  p^  is  one  of  the 
state  vector  initial  conditions,  then 


3yk 

*Pi 


= for  p.  i y£ 


t = t 


t = t 


In  general  a change  in  p^  will  affect  the  both  through  an  explicit 
dependence  of  F^  upon  p^,  and  also  implicitly  through  a change  evoked  in 
the  state  vector  y by  the  integrated  effect  of  the  p^  perturbation. 
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Namely, 

3F, 


3Pi  _ 
t,p^i 


3Fk  A 9Fk  3yji 

3—  + 22  * 3——  k = 1,3  (14) 

9Pi  - _ £-1  3y£  3Pi 

t,p  t i,y  1 t,p,y  j*£ 


The  vector  SF^/ 3p^ | t - -^.is  found  through  differentiation  of  the  explicit 
dependence  of  the  upon  p^,  therefore  when  p^  is  one  of  the  libTation 
initial  conditions  this  term  vanishes.  From  equations  (4)  we  have 


3F1 

3p. 

esc  0 | 

^sin  <{) 

3pT  + 

3F2  _ 
3p. 

cos  <f> 

3u>^ 

3p. 

sin  <p 

3F3  _ 

3(L3 

cos  0 

3F1 

3p. 

3?i  ’ 

3p. 

To  evaluate  the  quantities  Sto^/Sp^  we  may  substitute  equations  (9)  into 
equations  (2)  and  then  differentiate: 

*1  3a  V rJU  _3_/3Ub\  J_/?\YLt  3 MU 

35T  * "Vs  357  *^j9\  A [_x2  3p.  Ux3/  ‘ x3  3Pi  \3x2/J  Tbl  Jp“'  k)f 

157  - xi  *T“4r(r)} 

157 - -“i“2  ■§:  ‘ £ dl  T [X1  3?7  (3J7)  - x2  357  (3^7)]  *Tb3  357  (e)}  C16> 


The  partial  derivatives  with  respect  to  specific  parameters  in  equations  (16) 
are  simple  and  can  be  found  in  Appendix  A.1. 


■ - mu  - m irnmmmmMm* 


The  only  quantity  in  (14)  which  remains  to  be  formulated  is  the  3x6 
"feed-back"  matrix  SF^/ay^;  it  is  this  term  which  supplies  the  initial 
condition  partials.  First  we  differentiate  the  defining  equations  (4) 
for  the  F^  with  respect  to  each  component  of  the  state  vector: 


3F3  ^3  3F1 

— = — - cos  0 — 

dip  dip  dp 


(17) 


To  obtain  the  quantities  3u^/3y£  we  may  again  substitute  equations  (9) 
into  equations  (2)  and  differentiate,  this  time  with  respect  to  the 
Euler  angles: 
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f I 


9u)7 


B.  w 


J1  W 


♦£  % 

b=®,0  B 


+ (j 


3u  j 

3 Wo 


3 3y„ 


9ub 

-x  JL| 

3*1 30b  j 

\«V 

3xl 

1 9>'* 

V3V 

' 

3oj7 

WQ 


3uj 


3u>, 


'Y  “l  W ‘ “2  W. 


,y  \|  j 
b^»,8 c n 


3 f9ub 


3x.  3U, 
1 b 


SV\  3*2  3Ub 


\*qr 


Notice  that  the  summation  terms  are  zero  for  l - 4,  5,  or  6 since  the 
torques  are  not  dependent  upon  the  angular  velocity  of  the  Moon.  The 
terms  multiplying  the  moment  of  inertia  ratios  are  kinematic  and  can  be 
derived  from  equations  (3) : 


3a)j 

If 


3ui, 


= 0 


-jjg-  = sin  p cos  9 

3w. 

1 • . • 

= -0sm<j)  + p sin  9 cos  <(> 

3o)j 
3\p 

30 

3t0i 

3$ 


sin  4>  sin  0 


cos  p 


3(^2 

If 

3o>2 

~W 

3w2 

3 <p 

3oj2 

3^ 

9oj2 

30 

3^2 

dp 


= 0 


s \p  COS  <p  cos  0 


3o)3 

3^ 

3a)3 

30 

3oj, 


= 0 


-P  sin  9 


-0  cos  p - p sin  <j>  sin  0 = 0 


cos  <p  sin  9 


= -sin  <f> 


= 0 


dip 

3fa)3 

30 

3^3 

dp 


cos  9 


(19) 
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Inside  the  summation  we  have  the  terms  3x^/3y^  which  represent  the  change 
in  the  selenocentric  coordinates  of  a perturbing  body  with  respect  to 
changes  in  the  Euler  angles.  The  selenocentric  coordinates  are  trans- 
formed from  the  inertial  system  by  a rotation  matrix  R defined  by 


x = R 9,  4>)  X 


and  given  by  Goldstein  (1950)  as: 


ROf'.Q*^)  3 I"  cos  <p  cos  ip  - sin  4>  cos  0 sin  ij/  cos  4>  sin  ip  + sin  4>  cos  0 cos  ip  sin  4>  sin  8~ 

-sin  $ cos  ip  - cos  $ cos  0 sin  ip  -sin  $ sin  ij>  * cos  $ cos  0 cos  ip  cos  <p  sin  9 (20) 

. sin  0 sin  ip  -sin  0 cos  \ p cos  0 


Now  £ does  not  depend  on  the  orientation  of  the. selenocentric  coordinate 
system,  so : 


= (-cos  <J>  sin  ip  - sin  <j>  cos  0 cos  ip  ) £ + (cos  4>  cos  4>  - sin  <j»  cos  0 sin  i>  ) 


sin  <p  sin  0 sin  ip  - s in  sin  9 cos  ip  ^ + sin  4>  cos  9 £3 


(-sin  $ cos  ij; -cos  4>  cos  9 sin  <p  ) + (-sin  4>  sin  ij>  +cos  <j>  cos  6 cos  4> ) 


+ cos  sin  9 g, 


* (sin  4>  sin  ^ -cos  4>  cos  0 cos  4/ )£j  + (-sin  4>  cos  ij/ -cos  4>  cds  0 sin  40C-, 


cos  41  sin  9 sin  ip  Sj  - cos  <p  sin  9 cos  ip  C2  * cos  4>  cos  0 £3 
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a2x 

3xk  3xi 


(xl  6lk  * x2  S2k)(x2  {li  - X1  S2i>  4 (X1  5li  *x2  52i)Cx2  Slk  ' X1  42k5 


: 2 2.2 

(Xj  + x2) 


This  completes  the  derivation  of  the  variational  equations  for  our 
libration  model.  If  there  are  n parameters  for  which  we  desire  derivatives 
we  can  integrate  the  6n  first-order  variational  equations  in  parallel  with 
the  6 equations  of  motion  by  any  suitable  numerical  method. 


VERIFICATION 

The  equations  of  motion  and  variational  equations  have  been  integrated 
using  an  Adams-Moulton  predictor-corrector  method  (Smith,  1968)  at  8 steps/day 





mi 


- 
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within  the  framework  of  the  Nl.I.T.  Planetary  Ephemeris  Program  (Ash,  1965). 

The  orbital  motions  of  the  Moon  and  Earth  were  read  from  ephemeris  tapes 
which  were  integrated  from  initial  conditions  obtained  in  a least-squares 
fit  to  data:  5 years  of  laser  ranging  data  of  the  Moon  (King  et  al.,  1976), 

and  many  years  of  optical  and  radar  observations  of  the  planets  (Ash  et  al., 
1971).  The  variational  equations  have  proven  to  be  consistent  with  the 
equations  of  motion  to  better  than  .01%  by  finite  differencing  of  the 
latter;  the  fact  that  one  may  check  variational  equations  in  this  manner 
is  another  strong  incentive  (in  addition  to  the  operational  convenience  of 
the  parallel  integration  of  the  partial  derivatives)  for  their  use. 

In  order  to  verify  the  accuracy  of  our  integration  of  the  equations 
of  motion,  we  have  compared  the  integrated  Euler  angles  with  LLB-5.  The  LLB-5 
Cassini  angle  initial  conditions  were  transformed  to  the  corresponding 
Euler  angle  initial  conditions  referenced  to  the  inertial  1950.0  mean 
equator-equinox  system.  The  lunar  gravitational  harmonic  coefficients, 
moment  of  inertia  ratios,  and  the  mass  of  the  Earth  were  also  set  to  the 
values  used  to  generate  LLB-5.  The  result  of  the  integration  is  compared 
with  LLB-5  in  Fig.  2,  where  differences  in  the  Cassini  angles  (which  must 
be  transformed  back  from  the  Euler  angles  of  PEP's  output  at  each  tabular 
point)  are  plotted  as  functions  of  time.  The  greatest  discrepancy  is  in 
T,  and  has  the  form  of  a sinusoid  with  a 2.9  year  period  and  an  amplitude 
of  0.29”.  This  period  is  that- -of  the  free  mode  of  libration  in  longitude, 
corresponding  to  a homogeneous  solution  of  Euler's  equations.  This  mode 
has  been  stimulated  by  a longitude  offset  of  0.29”  between  the  lunar 
ephemeris  used  to  create  LLB-5  and  that  of  King  et  al.  The  cause  of  this 
orbital  offset  is  unknown,  but  under  investigation. 
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In  order  to  minimize  the  effects  of  the  JPL-MIT  mean  lunar  orbit 
difference  for  our  comparison  of  the  two  lunar  libration  integrations, 
we  introduced  three  ad  hoc  small  angles  to  describe  an  arbitrary  orienta- 
tion difference  between  coordinate  systems;  specifically,  we  modelled 
the  angles  as  constant  biases  in  each  of  the  three  Cassini  angles. 

Since  the  Cassini  angles  describe  the  difference  between  the  true  lunar 
libration  and  the  mean  lunar  orbit  (vis-a-vis  Cassini's  laws),  a bias 
in  a Cassini  angle  follows  directly  from  a corresponding  offset  in  the 
orientation  of  the  mean  lunar  orbit.  These  biases  were  simultaneously 
fit  with  six  initial  conditions  of  rotation  to  minimize  the  sum-squared 
differences  of  Cassini  angles  between  the  PEP  integration  and  LLB-S. 

The  adjusted  initial  conditions  were  then  integrated  and  the  new  Cassini 
angle  differences  from  LLB-5  are  plotted  with  biases  included  in  Figure  3; 

the  estimated  biases  and  post-fit  rms  Cassini  angle  differences  (about 
the  biases)  are  displayed  in  Table  1. 

The  largest  differences  between  the  PEP  post- fit  integration  and 

LLB-5  appear  about  equally  in  the  Cassini  angles  p and  la,  and  correspond 
to  a lunar  surface  displacement  of  about  16  cm  from  each  angle,  or  about 
4 cm  in  range  to  the  retro-reflectors  farthest  from  the  sub- Earth  point. 
The  sinusoidal  signatures  are  predominantly  those  of  the  free  libration 
modes  of  p and  la  — a combination  of  27.2  day  and  24  year  components 
in  both.  These  free  librations  are  of  a magnitude  consistent  with  the 
hypothesis  that  they  are  being  stimulated  by  the  known  differences  in 
the  lunar  ephemerides  used  in  the  integrations. 

In  order  to  ascertain  which  libration  model  more  closely  represents 
the  true  motion  of  points  on  the  lunar  surface  about  the  center  of  mass 
we  have  performed  parallel  fits  to  laser  ranging  observations  of  the 


Apollo  retro-reflectors.  The  maximum  difference  between  the  Euler  model 
solution  and  a corresponding  solution  using  LLB-5  was  0.2  ns  (3  cm  in 
range)  rms  for  a series  of  range  observations  on  a specific  retro-reflector 
over  a period  of  2 or  3 years,  with  the  rms  residual  for  the  Euler  model 
equal  or  lower  than  LLB-S  for  all  such  series.  It  should  be  emphasized, 
however,  that  the  lunar  orbital  emphemeris  employed  in  these  data  analyses 
was  the  same  ephemeris  which  was  used  to  generate  the  Euler  angle  libration 
model  and  perhaps  the  lower  residuals  merely  reflect  the  consistency  of 
orbital  and  rotational  models. 

In  our  development  of  the  Euler  angle  model  we  have  made  many  approx- 
imations, some  of  which  we  will  continue  to  ignore  and  others  we  will 
correct  in  the  future.  The  lunar  orbital  and  rotational  motions  are  cross- 
coupled  and  we  are  in  the  process  of  modifying  our  computer  algorithms  to 
integrate  them  simultaneously.  The  Earth’s  oblateness  appears  to  affect 
libration  at  above  the  .01"  level,  so  we  plan  to  include  the  low-order 
gravity  field  of  the  Earth  in  our  model.  Finally,  the  Moon  is  not  a rigid 
body,  and  is  subject  to  solid  body  tides  raised  mainly  by  the  Earth, 
resulting  in  an  inertia  tensor  that  is  non-diagonal  and  time-variable;  also 
there  is  undoubtably  some  tidal  dissipation  in  the  Moon.  The  magnitudes 
of  these  effects  are  difficult  to  estimate,  depending  as  they  do  upon 
unknown  quantities  such  as  the  lunar  Love  numbers,  and  Q for  the  Moon, 
though  a rough  calculation  places  their  effect  at  about  the  level  of  the 
other  errors  in  our  model,  so  we  plan  to  incorporate  a flabby  Moon  into  our 
theory  as  well.  Ultimately,  we  hope  to  reduce  the  total  libration  errors 
to  below  the  anticipated  future  lunar  laser  range  uncertainties  of  3 cm, 
so  that  the  full  benefit  of  the  quality  of  the  laser  ranging  data  can  be 


realized. 
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APPENDIX  A.  1 

PARTIAL  DERIVATIVES  WITH  RESPECT  TO 
SPECIFIC  LIBRATION  PARAMETERS 

In  order  to  derive  the  partial  derivatives  indicated  in  equation  (16) 
we  must  first  determine  which  members  to  include  in  our  set  of  adjustable 
libration  parameters.  In  principle  any  unknown  parameter  affecting  the 
libration  should  be  modelled  as  such  and  have  derivatives  generated  for 
it;  practically,  we  need  only  model  the  small  set  of  parameters  whose 
uncertainties  influence  the  libration  to  a measurable  extent  (defined, 
for  the  purpose  of  this  paper,  as  greater  than  1 part  in  10^  of  the  whole 
libration).  We  are  left  with  the  following  (significant)  parameters. 

1)  six  initial  conditions  of  the  state  vector 

2)  J2 

3)  £5  and  y 

4)  third  and  higher  order  harmonic  coefficients 

Note  that  we  have  arbitrarily  chosen  J2,  B,  and  y as  our  independent 
second-degree  coefficients.  C22  is  a combination  of  all  three,  while 
a depends  only  on  0 and  y,  and  since  we  choose  the  lunar  principal  axes 
of  inertia  to  define  our  selenocentric  coordinate  system,  C2^,  S2^,  and 
S22  are  identically  zero. 

Now  we  will  derive  the  derivative  terms  in  equation  (16).  38/3p^ 

and  Dy/Spi  are  non- zero  only  when  pi  = 0 or  y respectively.  3a/3pit 
SA’Vsp^,  3B"V3p^,  and  3C  */3p^  are  all  zero  for  p^  other  than  the 
second  degree  parameters  0,  y,  and  J^,  when  we  get  the  following: 
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The  remaining  chain-rule  terms  in  equations  (16)  are  the  3/3p  (3U  /3x,  ) 

i b k' 

The  harmonic  coefficient  partials  are  merely  the  terms  of  the  linear 


combination  which  the  coefficients  multiply: 
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For  pA  one  of  0,  y,  or  J we  apply  the  chain  rule  to  find  the  (additional, 
in  the  case  of  J2)  change  in  the  potential  gradient  due  to  the  change 
induced  in  C^: 
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The  second  term  in  the  RHS  of  equation  (A-3)  is  given  by  (A-2b)  with 
n * m = 2. 
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Comparison  of  Best-fit  PEP  Euler  Integration  with  LLB-5 


bias 

post-fit  rms  difference 
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